In [ ]:
import io
import math
import os
import pathlib

import holoviews as hv
import hvplot.pandas
import numpy as np
import pandas as pd
import pymap3d
import xarray as xr

hv.extension("bokeh")
np.set_printoptions(suppress=True)
No description has been provided for this image No description has been provided for this image
In [ ]:
def get_skews_and_base_cfls(lons, lats, depths) -> np.ndarray:
    # The shape of each one of the input arrays needs to be (3, <no_triangles>)
    #ell = pymap3d.Ellipsoid.from_name("wgs84")
    ell = pymap3d.Ellipsoid(6378206.4, 6378206.4, "schism", "schism")
    local_x, local_y, _ = pymap3d.geodetic2enu(lats, lons, depths, lats[0], lons[0], depths[0], ell=ell)
    areas = (local_x[1] * local_y[2] - local_x[2] * local_y[1]) * 0.5
    rhos = np.sqrt(areas / np.pi)
    max_sides = np.maximum(
        np.sqrt(local_x[1] ** 2 + local_y[1] ** 2),
        np.sqrt(local_x[2] ** 2 + local_y[2] ** 2),
        np.sqrt((local_x[2] - local_x[1]) ** 2 + (local_y[2] - local_y[1]) ** 2),
    )
    skews = max_sides / rhos
    base_cfls = np.sqrt(9.81 * np.maximum(0.1, depths.mean(axis=0))) / rhos / 2
    return skews, base_cfls

def get_skews_and_base_cfls_from_path(path: os.PathLike[str] | str) -> np.ndarray:
    ds = xr.open_dataset(path, engine='selafin')
    tri = ds.attrs['ikle2'] - 1
    lons = ds.x.values[tri].T
    lats = ds.y.values[tri].T
    depths = - ds.B.isel(time=0).values[tri].T
    skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)
    return skews, base_cfls
In [ ]:
file = "/home/tomsail/Documents/work/models/meshes/slf/v1p2.slf"
file = "/home/tomsail/Documents/work/models/meshes/slf/v2p1.slf"
ds = xr.open_dataset(file, engine='selafin')
skews, base_cfls = get_skews_and_base_cfls_from_path(file)
In [ ]:
CFL_THRESHOLD = 0.4
for dt in (1, 50, 75, 100, 120, 150, 200, 300, 400, 600, 900, 1200, 1800, 3600):
    violations = (base_cfls * dt < CFL_THRESHOLD).sum()
    print(f"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%")

   1      5887124   100.00%
  50      1566579    26.61%
  75      1095580    18.61%
 100       824343    14.00%
 120       664172    11.28%
 150       344840     5.86%
 200       158392     2.69%
 300        61959     1.05%
 400        11172     0.19%
 600         3120     0.05%
 900         1934     0.03%
1200          522     0.01%
1800           11     0.00%
3600            0     0.00%
In [ ]:
pd.DataFrame({"skew": skews}).describe()
Out[ ]:
skew
count 5.887124e+06
mean 2.826159e+00
std 1.493293e-01
min 2.506716e+00
25% 2.724916e+00
50% 2.790498e+00
75% 2.894929e+00
max 2.021513e+01
In [ ]:
df = pd.DataFrame({"cfl": base_cfls * 400})
df[df.cfl < 0.4].describe()
df[df.cfl < 0.4].hvplot.hist()
Out[ ]:
In [ ]:
tri = ds.attrs['ikle2'] - 1
nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=["lon", "lat", "depth"])
elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=["no_nodes", "n1", "n2", "n3"])
elements = elements.assign(base_cfl=base_cfls)
elements.head()
Out[ ]:
no_nodes n1 n2 n3 base_cfl
0 3.0 802081.0 65536.0 2267476.0 0.002489
1 3.0 585525.0 65536.0 802081.0 0.002548
2 3.0 65536.0 585525.0 780858.0 0.002697
3 3.0 251963.0 65536.0 780858.0 0.002613
4 3.0 65536.0 251963.0 820904.0 0.002469
In [ ]:
min_cfl_per_node = pd.concat([
    elements[["n1", "base_cfl"]].groupby(["n1"]).base_cfl.min(),
    elements[["n2", "base_cfl"]].groupby(["n2"]).base_cfl.min(),
    elements[["n3", "base_cfl"]].groupby(["n3"]).base_cfl.min(),
], axis=1).min(axis=1)
min_cfl_per_node.head()
Out[ ]:
1.0    0.005506
2.0    0.007554
3.0    0.003865
4.0    0.018533
5.0    0.021136
dtype: float32
In [ ]:
dt = 200
df = nodes.assign(
    cfl=min_cfl_per_node * dt,
    # CFL_violation nodes have a value of 1 if there is no violation and 4 if there is a violation. 
    # We do this in order to plot the points with a different size
    cfl_violation=((min_cfl_per_node * dt < CFL_THRESHOLD) * 3) + 1   
)
df.head()
Out[ ]:
lon lat depth cfl cfl_violation
0 -180.000000 90.000000 -4228.0 NaN NaN
1 134.560318 -65.732597 -107.0 1.101187 1.0
2 -13.422298 29.300503 -68.0 1.510851 1.0
3 12.516141 54.809700 -20.0 0.772975 1.0
4 131.169098 -0.395918 -300.0 3.706690 1.0
In [ ]:
plot = df[df.cfl_violation == 4].hvplot.points(
    'lon', 
    'lat',
    c="depth",
    cmap="jet",
    geo=True,
    tiles="EsriImagery",
).options(
    width=1200, height=900
)
len(df[df.cfl_violation == 4])
Out[ ]:
135053
In [ ]:
plot
Out[ ]: